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The foundation for an accurate and unifying Fourier-based theory of diffusion weighted 
magnetic resonance imaging (DW-MRI) is constructed by carefully re-examining the first 
principles of DW-MRI signal formation and deriving its mathematical model from scratch. 
The derivations are specifically obtained for DW-MRI signal by including all of its elements 
(e.g., imaging gradients) using complex values. Particle methods are utilized in contrast to 
conventional partial differential equations approach. The signal is shown to be the Fourier 
transform of the joint distribution of number of the magnetic moments (at a given location 
at the initial time) and magnetic moment displacement integrals. In effect, the /c-space is 
augmented by three more dimensions, corresponding to the frequency variables dual to 
displacement integral vectors. The joint distribution function is recovered by applying the 
Fourier transform to the complete high-dimensional data set. In the process, to obtain a 
physically meaningful real valued distribution function, phase corrections are applied for 
the re-establishment of Hermitian symmetry in the signal. Consequently, the method is 
fully unconstrained and directly presents the distribution of displacement integrals without 
any assumptions such as symmetry or Markovian property. The joint distribution function 
is visualized with isosurfaces, which describe the displacement integrals, overlaid on the 
distribution map of the number of magnetic moments with low mobility. The model 
provides an accurate description of the molecular motion measurements via DW-MRI. The 
improvement of the characterization of tissue microstructure leads to a better localization, 
detection and assessment of biological properties such as white matter integrity. The 
results are demonstrated on the experimental data obtained from an ex vivo baboon brain. 
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1. INTRODUCTION 

Since the conception of mathematical models for the effect of 
the magnetic moment diffusion in nuclear magnetic resonance 
(NMR) experiments by Hahn (1950), Carr and Purcell (1954), 
and Torrey (1956), several methods have been proposed for 
analysis of diffusion-weighted (DW) magnetic resonance imag- 
ing (MRI) signal. These advancements endowed with the non- 
invasive, in vivo nature of the technique, have propelled the initial 
utilization of DW imaging measures, e.g., apparent diffusion coef- 
ficient in early detection of ischemia (Moseley et al., 1990; Baird 
and Warach, 1998), to many highly crucial areas in research and 
clinical imaging: for example in cancer diagnosis (Song et al., 
2002; Turkbey et al, 2009; Xu et al, 2009), follow-up on treat- 
ment, pre- and post-operative assessment for different organs 
[e.g., fiber tracking (Conturo et al., 1999; Mori and van Zijl, 
2002)] white matter integrity assessment (Budde et al., 2007; 
Correia et al, 2008) as in monitoring of neurological diseases such 
as multiple sclerosis (Song et al., 2002) and disorders (Ciccarelli 
et al, 2008) like schizophrenia (Seal et al., 2008; Voineskos et al., 
2010) and Alzheimer's disease (Mielke et al., 2009), as well as 
neonatal development (McKinstry et al., 2002) and traumatic 
brain injury (Mac Donald et al., 2011). 



In brief, diffusion weighted magnetic resonance imaging 
(DW-MRI) has become an indispensable and versatile technique 
playing an important role in several applications by its ability to 
estimate diffusion. The abundance of DW-MRI models is an indi- 
cator of room for improvement as well as the necessity for unifi- 
cation [see Ozcan et al. (2012) for a detailed account of the partial 
differential equation (PDE) based adaptation's implications as 
well as a thorough mathematical analysis and a description of the 
background of existing methods] . 

DW-MRI's aim is to obtain measures and characterization 
of microstructure by investigating the diffusion process. Several 
methods and models have been proposed, all originating from 
the seminal work of Stejskal and Tanner (1965). Therein, under 
the influence of the additional motion sensitizing magnetic field 
gradients, the self-diffusion PDE of the magnetic moments is 
included in the Bloch PDE to model the attenuation in the 
DW-NMR spectroscopy signal. The result is the estimation of 
the scalar diffusion coefficient of the entire sample. In a sense, 
DW-NMR added another dimension, i.e., the magnetic moment 
motion, to the spectroscopic information even before the intro- 
duction of magnetic moment position later by the invention of 
MR imaging. 
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Accordingly, herein, DW-MRI model is naturally progressed 
to a higher dimensional construct that jointly presents magnetic 
moment position and motion. This is achieved by carefully re- 
examining the first principles of DW-MRI signal formation using 
particle methods in the spirit of the work of McCall et al. (1963). 
The mathematical model constructed in section 2.1 is specifically 
obtained for DW-MRI signal (rather than DW-NMR) by includ- 
ing all of its elements (e.g., imaging gradients) using complex 
values without taking the signal's magnitude. 

The approach reveals that for an / mr -dimensional MRI slice, 
the DW-MRI complex valued signal that comes out of the scan- 
ner is the (Z mr + 3) -dimensional Fourier transform of the joint 
distribution function of the number of magnetic moments (that 
are at a given position at the initial time) and their displacement 
integrals. In other words, the first Z mr dimensions correspond 
to the usual MRI /c-space with position information and the 
remaining three dimensions constitute the frequency space of 
displacement integrals. The values of imaging and motion sen- 
sitizing magnetic field gradient vectors together define in the 
(l mr + 3) -dimensional Fourier space the sampling points of the 
joint distribution function's Fourier transform. The distribution 
function is recovered by taking the Fourier transform of the com- 
plete data directly (i.e., without any scaling or use of magnitudes), 
giving the method its name: Complete Fourier Direct (CFD) 
MRI (Ozcan, 2010b). 

2. MATERIALS AND METHODS 

2.1. COMPLETE FOURIER DIRECT MRI SIGNAL FORMATION 

The MR signal is generated by the vectorial sum of transverse 
magnetization of magnetic moments (Haacke et al., 1999): 

M(t) = Y J m{t). (i) 

i 

By neglecting the effect of spin-spin relaxation, the evolution 
of the transverse magnetization of the ith magnetic moment 
is described in a standard manner by a rotating magnetization 
vector according to Bloch equations (see Appendix B): 

ttii(t) = exp(-j y £2,-) ff!,-(fo). (2) 



Here, y is the gyromagnetic ratio, the transverse magnetiza- 
tion vector, m,-, is written in complex number form with m,(to) 
denoting the initial magnetization tipped to the transverse plane, 

f2, = j G(xi, x) ■ xi(x)dx (3) 

to 

describes the phase (when multiplied by y) as a function of the 
magnetic field gradients G(x, t) G R 3 , and the position of the 
magnetic moment is x; e R 3 . The time-dependent position, x,-, 
in Equation (3) affects the phase, f2,-, thereby also affecting the 
total signal in Equation (1). 

Any kind of displacement (such as Brownian motion, molec- 
ular movement in biological tissue with different medium and 
obstacles, coherent motion or any combination thereof) is incor- 
porated straight into the model, by modeling the position in a 
general and direct form herein without any stochastic assump- 
tions [such as Markovian property used in Wedeen et al. (2005)] 
on the motion 

Xi(t) =Xi(t 0 )+Wi(t), (4) 

where w,-(f) e K 3 represents the displacement of the magnetic 
moment from its initial position, x;(fo), [i.e., w;(fo) = 0]. The 
only physical requirement is the continuity of w;(f) since a mag- 
netic moment cannot disappear at a given point and reappear at 
another. 

The signal is calculated using Equations (1- 4) in reverse 
order. Following the motion described by Equation (4), the phase 
of the ith magnetic moment in Equation (3) during the dig- 
ital acquisition period of the two dimensional imaging (Z mr = 
2) pulsed-gradient spin-echo (PGSE) sequence of Figure 1, is 
obtained after tedious but routine derivations [see Appendix B 
for a brief exposition of the derivations and Ozcan (2012)] using 
the definitions of the variables in Figure 1 with G* € R 3 denot- 
ing the magnetic field gradient vectors labeled as read-out, ro; 
phase encode, pe; slice select, ss; and diffusion, D. Omitting rou- 
tine calculations for trapezoidal shapes for clarity, the derivation 
is carried out assuming ideal gradient amplifiers providing rect- 
angular shaped gradient pulses. The initial time, to, is the end time 
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FIGURE 1 | The pulsed-gradient spin-echo (PGSE) pulse sequence and 
the definition of the variables used in the calculations. SS-EX is for the 
slice select gradient during the excitation (n/2) pulse, RO for read out, PE for 
phase encode, SS for is the slice select gradient, Diff marks the motion 



sensitizing pulses, SS-PI is the slice select gradient during the u-pulse 
and ACQ stands for digital acquisition period. In practice, the MR pulse 
sequences implement the rewind (rw) gradients such that the amplifiers 
are turned on and off at the same times. 
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of 7t/2 radio frequency (RF) pulse when the magnetization is fully 
tipped to the transversal plane. The resulting phase of the trans- 
verse magnetization is a function of time, f , and, the imaging and 
diffusion gradients (see Appendix B): 

&i(t, Gpe, G ss , Gd) 

= ((f - facq - Af rw ) G ro - Af rw Gpe) • Xj(f 0 ) - Af rw G ss • X,(t 0 ) 

(5) 

+G D • W d + ((t d4 - f d3 ) - (tdi - «ai))G D ■ Xi(to) + 4>«i (6) 



+G ro • wf q (t) - (G ro + G pe + G ss ) • W? 



(7) 



The second term in Equation (6) removes the injection of the 
initial position into the DW signal because of equal pulse dura- 
tion times, 8 = f^4 — tfo = — ftfi- The term tt,,; describes 
the systematic phase (see Appendix Bl) created by the it-pulse 
slice select gradient (SS-PI) and in this work it will be auto- 
matically taken out by the phase correction algorithm in sec- 
tion 2.2. Equations (6) and (7) incorporate three integrals of the 
displacement w,(f) : (W d , W acq , W™) corresponding to the dis- 
placement integrals for diffusion (d), analog to digital conversion 
acquisition (acq), and initial rewind (rw) gradient time periods, 
respectively 



j wi(x) dx — j 1 



Wi(x) dx, 



[ irvi 

Wf q (f) = J w,(x)dx,W™ = J w,(x)dx. 



(8) 



(9) 



First term in Equation (5) is the definition of the fc-space in 
regular MRI, fc mr e R 2 : 



affine function of fc n 



[~gn 



2 /Afrw, 



3] 



(12) 



Consequently, the phase, fi,*, of Equation (3) is expressed con- 
cisely by defining kv = Gd : 



£2; = k mr ■ x,(t 0 ) + k D -W° + k r 

■icq 
1,1 



WI 



+ W- cq (k mi i)g m i + tp slke 



(13) 



reflecting the effect of initial position and displacement integrals 
on the phase 2 of each magnetic moment. Since <p s ii ce is constant 
for all i, it is taken out of Equation (13) with the appropriate 
rotation of the magnetization coordinate frame on a slice by slice 
basis. 

Finally, assuming that all of the magnetic moments have 
the real valued initial magnetization m,(fo) = mo, Equation (1) 
can be re-written using Equation (13) to reveal a Fourier 
relationship, 

S c fd(kmr, fc D, krw) = m 0 ^ exp(-; Y ^i(fcmr, fc>, W)- (14) 



A more efficient way to evaluate the sum in Equation (14) is 
first to group the magnetic moments with the same position- 
displacement properties and then to count the numbers elements 
in the groups. 

Definition 1. The joint position-displacement integral distribu- 
tion function, P l0 ^ al (jc, W), is defined as the number of mag- 



the initial position isl at 
possessing the displacement integral values of W = (W 



netic moments with 



wf cq ) € : 



time to, 

d jyrw 



fcmrO, Gpe) = (f - t, 



acq 



Af rw ) G ro — Af rw G 



pe 



with the additional term, 



<Pslice 



-Af rw G ss • Xi(t 0 ), 



(10) 



(11) 



which is constant because the slice select axis component of X;(frj) 
is the slice position 1 . 

Without loss of generality, by adopting the imaging coor- 
dinate frame defined by the directions of the read-out, phase 
encode and slice select gradients, G ro = [g ro i, 0, 0], G pe = 
[0, g pe 2, 0], G ss = [0, 0, gssi], time andg pe 2 become functions 
of fcmr = [fcmrl, fcmr2, k m[3 ] using Equation ( 10): 

t = fcmrl/grol + facq + Af rw and g pe2 = -fcmrt/Afrw. 

Accordingly, W acq (f) becomes a function of fc mr i and the coeffi- 
cients of W™ in Equation (7) are written as a vector which is an 



'The refocusing lobe of the slice select gradient after the RF pulse denoted 
by SS in Figure 1 adjusts ip sl i cc to be constant throughout the slice in the slice 
select direction. 



The signal in Equation (14) is calculated by integrating over 
the whole position-displacement space (absorbing mo into PHj 
for ease of notation): 

Scfd(fcmr, to, fcrw) = f Pgffc W) 

x exp(-7 Y (fcmr ' X + fc D • W d + fc rw • W rw 
+ Wl cq {k mrl )g ml ))dxdW (15) 

= ^{P^ al }(fc mr ,fc D ,fcrw). (16) 

Equation ( 15) is by definition the Fourier transform of P^jr with 
non-linearities added by W* cq . 

Among the elements of W, the focus is on the most descriptive 
MRI observable, W d . Its distribution is obtained by marginalizing 



"Jrr., e , . magnetic field X time i . r i j r • magnetic field 

The unlt for fc mr is dlstance but for k D and k IVI it is g istance . 

They are both in accordance with the units of the position and the displace- 
ment integrals in Equation (13). After k c [& is multiplied by the gyromagnetic 
ratio y, with unit (magneticfield x time) -1 , the product becomes unitless in 
Equation (14). 
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W rw from Equation (15) 

P cfd (x, W d ) = f Prffix, W d , W™)dW™ 

= ^i.to)f S cfd(/c m r,/c D ,0)}. (17) 

However, the affine dependence of /c rw in Equation (12) makes it 
impossible to fix fc rw = 0 and to sample in (/c mr , /cd, 0) subspace. 
The following example demonstrates how the affine dependence 
affects the measurements by using a two dimensional Gaussian 
function, exp (— (k\ + k\)), with the "undesirable" variable k 2 
sampled on a line k2 = a k\ + b : 

^{^(-(kl + iah + b) 2 ))} 



2^71(1 + fl 2 ) 



exp 



<j + J Aabxi + 4b 2 
4a 2 + 4 



The result is complex valued in comparison to the real valued 
Fourier transform of exp(— k\) which can be obtained by setting 
a = b = 0. 

2.2. PHASE CORRECTIONS FOR THE ESTIMATION OF P cfd 

In addition to inherent affine dependence and non-linearities, 
different experimental factors, noise, hardware imperfections 
etc., affect the DW-MR signal adversely. CFD-MRI addresses 
these issues by adopting a pivotal, physically meaningful 
standpoint originating from the following Fourier transform 
property (Bracewell, 2000): 

(real valued function) t> J- t> (Hermitian symmetric function). 

Accordingly CFD-MRI reconstruction is based on the 
following: 



Since by definition P^ 3 ' is real valued, S c fd is Hermitian symmetric. 



Furthermore, an immediate implication of the transform 
property and the Hermitian symmetry of S c fd is that theoret- 
ically, taking its magnitude before Fourier transforming will 
result in a symmetric real valued distribution under ideal condi- 
tions. As noted above, in practice the experimental S c fy is never 
Hermitian symmetric resulting in an asymmetric magnitude. 
Consequently, the magnitude's Fourier transform used in exist- 
ing methods (Callaghan, 1991; Wedeen et al., 2005), results in 
a complex valued (Hermitian symmetric) distribution function. 
The difficulty of a physical interpretation forced those methods 
to take the magnitude of the transform as well to obtain a real 
valued function. 

Herein, in order to obtain a real valued P^a, the re- 
establishment of Hermitian symmetry in the signal during the 
computation of the inverse transform for Equation (17), is real- 
ized by phase corrections. The strategy is similar in principle to 
the correction of the A: ml --space center's (echo time) shift during 
the read-out period of acquisition in MR imaging. The resulting 



linear phase shift in the physical read-out axis uniformly and 
systematically appears in all of the data. The shifts in both phase- 
encode (e.g., due to sample shaking) and read-out directions are 
corrected by first determining from the Fourier transform in k ml , 

£° mPleX ((Xro, *pe)> fc D) = FZ l {S c fd(fcmr, k D, Wl- (18) 

the angle (^/™ m P lex ) f rom the signal regions along the center lines 
of each physical direction at fco = 0, 

r ro (x ro ) « Z7™ mplex ((x ro , 0), 0), r pe (x pe ) « Z/" mplex ((0, x pe ), 0). 

(19) 

The phase corrections are then applied systematically at each 
value of /cd (see Figure 3): 



h m (x, k D ) = /™ mpleX (x, fc D ) exp(-J (rro(*ro) + f-pe(Xpe)))- 



The Fourier transform in the remaining variables, 
P cfd (x,W d ) = ^ D l {I kml (x, fc D )}, 



(20) 



(21) 



is evaluated sequentially in each /co-dimension with the aim of re- 
establishing the Hermitian property, h^ix, —ko) = 1% (x, k^), 
using the following steps. 

CFD Phase correction algorithm: 

1. Pick a pixel at location x c , preferably near the center of the 
image where tissue or a good signal area is located. 

2. Starting from the first direction, I = 1 of the fco space cal- 
culate the phase on the line passing through the origin 
(i.e., [fc D i, 0, 0], [0, k D2 , 0], [0, 0, A: D3 ], respectively for 1 = 
1,2, 3), e.g., /.I km (x c ,(k m ,0,0)). 

3. Investigate the plot of the phase versus fcrj;. Pick as many 
as possible consecutive values of k^i near 0 without sudden 
changes to assure high signal to noise value. 

4. Construct a polynomial of degree m (with m less than the 
number of points) that approximates the phase at the selected 
points. The polynomial's constant term must be set to be 0 
to guarantee thatlj; mr (x c , 0) remains unchanged. For example, 
for the first direction, at selected values of fcoi the polynomial 
looks like 



lh(x c , (k m , 0, 0)) sa r m {k m ) = a m (k m f 



+ a m - 1 (k m ) m 1 + ... + aik 



Dl 



(22) 



as demonstrated in Figure 2. 

5. Apply the same phase correction systematically to the 
entirety of the data along the /th direction at each of the 
other dimensions, at all of the pixel locations. For exam- 
ple in the first direction, fcoi, update to be equal 
to h m (x, (fc D i, k m , km)) exp(-j r D1 (fc D i)) for all k D2 , k m 
and x. 

6. Repeat steps 2-5 for the remaining directions. 

The algorithm transfers the signal to the real channel by 
preserving its energy as the phase corrected spin-echo image 
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FIGURE 2 | Ozcan (2012): Top row: The Nyquist plots of uncorrected (left) 
and corrected (right) S^d obtained from the experimental data described 
in section 2.4. The plots show data acquired at each diffusion gradient value 
k D i on the complex plane. Bottom row: The magnitude and phase plots of 
the data. Uncorrected data (bottom row left, second column) exhibit a linear 



phase shift around 0 frequency, indicative of coherent motion. After the phase 
corrections obtained using the polynomial 0.266 fcci estimated from the points 
fc D i = —6, 0, 6, 12 Gauss/cm, the magnitude is unchanged but the signal's 
imaginary part is smaller for the corrected values visible by the difference 
between the vertical axis spans of Nyquist plots and the phase plots. 



without diffusion gradients, 1^ (x, 0) . The distribution of the 
magnetic moments with low mobility in all three directions [i.e., 
Pdd(x, 0)] shows the result of the transfer in Figure 3 for the 
sample described in section 2.4. 

In P c f,j (x, 0) , areas with high level of organization inducing low 
mobility, such as the corpus callosum (CC), the external capsule 
(EC), the mid-brain and the pons, appear brighter. The image is 
not an anisotropy map, e.g., mineral oil would appear brighter 
than water due to a smaller diffusion coefficient despite both liq- 
uids being isotropic. Spin-echo image is more blurred because it 
is a low pass filtered version of P^: 



7™ mplex (x, ho) = Fvl{Scfd(k mr , *d, fcrw)} =► h mi (x, 0) 



2 jr. f 

7 J 



Ptff\x, W) dW 



(23) 



(see Appendix C). 

CFD phase correction algorithm outperformed the fitting of 
the phase values up to the fourth degree multinomials in R 3 . 
The reasons behind this outcome, which will provide information 
about DW-MR signal artifacts, as well as inclusion of different 
functions for corrections will be investigated in the future. 

2.3. CFD-MRI SAMPLING AND WINDOWING 

Whereas the standard MRI field of view (FOV) calcula- 
tions (Haacke et al, 1999) are used for fc mr -space, the infinite 
bandwidth in ^D-space due to P c fd's finite support in W d -space 
(originating from finite length displacements) falls beyond the 



reach of the gradient hardware's limits for small diffusion gra- 
dient duration and separation times (8 and A, respectively in 
Figure 1). Even with a powerful gradient system, a large mag- 
nitude of /cd causes substantial signal uncertainties due to an 
increasing performance deterioration as the power requirements 
push the hardware to its limitations. 

With such a hardware constraint, in order to reduce ripple 
effects caused by truncation, P c fd's bandwidth (i.e., S c fd's sup- 
port) is shrunk by increasing 8 and A causing the dispersion 
(covariance) of the displacement integral W d (and therefore P c fd's 
support) to increase. This is directly visible in the special case 
of Brownian Motion characterized by the diffusion tensor D in 
Figure 4. P c fd and S c fd are zero mean Gaussians with covariances, 
respectively equal to (see Appendix A) 

E[W d (W d ) T ] = b t D and Q> t Dy l where b t = S 2 (A - 8/3) 

(24) 

(see Ozcan, 2009, 2010a) because the Fourier transform of a 
Gaussian with a covariance matrix D is also a Gaussian with 
covariance D~ 1 : 



jexp((w d ) 



D~ 



exp 



(klbk D ) 



(25) 



The procedure is graphically displayed in Figure 4 also empha- 
sizing the effect of ripples on the small values of P c fd which are 
especially important in revealing microstructure as explained in 
section 3.1. 
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FIGURE 3 | On the left, the fixed baboon brain images acquired in the 
anatomical-coronal plane. On the top row, the real and imaginary parts of 
/<r mr (x, 0) and on the bottom row, P C fd(x, 0) are displayed. The imaginary part 
is approximately 10% of the real part in both cases. On the right, full 



representation of P CTt j with isosurfaces (P c fd = 0.17, see section 3.1) around 
CC and EC junction. Starting from left bottom going clockwise, the sample 
pixels are from cerebro-spinal fluid (CSF), CC, white matter (WM) and CC, 
and EC junction, respectively (see also Figure 5). 



The second sampling criterion is an appropriate sampling rate 
i.e., sufficient number of points in fco-space to prevent aliasing 
artifacts on P c fy. This is constrained by the time available for 
acquisitions as each point in kj) -space requires the scan time of 
the entire /c mr -space. 

2.4. EXPERIMENTAL SETUP AND ANALYSIS METHODS 

A fixed baboon brain immersed in 4% paraformaldehyde 
was used for the experiments. The primate was prematurely 
delivered on the 125th day and sacrificed on the 59th day 
after delivery. All animal husbandry, handling, and procedures 
were performed at the Southwest Foundation for Biomedical 
Research, San Antonio, Texas. Animal handling and ethics were 
approved to conform to American Association for Accreditation 
of Laboratory Animal Care (AAALAC) guidelines. Further 
details of the preparation are described in Kroenke et al. 
(2005). 

The experiments were carried out on a 4.7Tesla MR scan- 
ner (Varian NMR Systems, Palo Alto, CA, USA) with a 15 cm 
inner diameter gradient system, 45 Gauss/cm maximum gradi- 
ent strength and 0.2 ms rise time using a cylindrical quadrature 
birdcage coil (Varian NMR Systems, Palo Alto, CA, USA) with 
63 mm inner diameter. CFD-MRI data were obtained using the 
standard pulsed-gradient spin-echo multi-slice sequence. The 
fc mr -space was sampled to result in images of 128 x 128 pixels 
with a FOV 64 x 64 mm 2 and 0.5 mm slice thickness. The /cd- 
space was sampled in a uniformly spaced Cartesian grid in a 
cube [—30, 30 Gauss/cm] 3 with 6 Gauss/cm sampling intervals 
at each dimension resulting in llxllxll voxels. The repeti- 
tion time Tr = 1 s, echo time Te = 56.5 ms, diffusion pulse sep- 
aration time A = 30 ms and diffusion pulse duration 8 = 15 ms 
were used. 



The data were transferred to a two quad core 2.3 GHz Intel 
Xeon® cpu and 8 GB memory Dell Precision Workstation 490 
running Windows XP® 64-bit operating system. The DWI data 
were placed in a 5-dimensional array in the computer mem- 
ory and the discrete Fourier transform (DFT) was computed 
along with the phase corrections. In-house Matlab® (Mathworks, 
Natick, MA, USA) programs were used for all of the computations 
and to display the graphics and maps. 

3. RESULTS 

3.1. VISUALIZATION OF THE CFD DISTRIBUTION 

The joint distribution's high-dimensionality [e.g., two dimen- 
sions for position in regular MR images (Z mr = 2), plus three 
more for displacement integrals] creates a visualization challenge 
which is addressed herein by using P c a(x, 0) as the background 
image. Furthermore, the isosurface 3 of normalized P c fd, 

P c a(x, W ) = 

Pcfd(x,0) 

is overlayed on the pixel at location x, as in Figure 3. For the 
sake of an objective assessment, the isosurfaces are defined using 
a common level value c (0 < c < 1), 

(W d 6l 3 : P cfd (x, W d ) = c} , 



3 Another approach in the literature is to present the value of the function 
over a sphere (Tuch, 2004; Wedeen et al., 2005). However, uniqueness is lost 
in this representation as demonstrated by these functions: f(x, y, z) = x 2 + 
y 2 + z 2 , h(x, y, z) = (f(x, y, z)) 2 , which have the same value on the unit 
sphere but different isosurfaces. 
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FIGURE 4 | Effects of fitting the bandwidth within the field of view (FOV) 
in the Fourier space on the left, with the covariance of the Gaussians 
given as D = (8 2 (A — 8/3) D) _1 . The sampling is done in the fixed interval, 
[—5, 5], of the Fourier Space followed by reconstruction on the right using the 



discrete Fourier transform. Theoretically this is equivalent to convolution with 
the sine function (see Brigham, 1988) in physical space. The ripple effects 
created by sine lobes diminish on the left as the Gaussian falls into the FOV 
with increasing S, A but constant D. 



over all locations. The key point is the choice of an appropriate 
c-value that will reveal the outskirts of P c fd corresponding to the 
small number of "scout" magnetic moments that travel further 
away thereby portraying the microstructure of the environment. 
In summary, 

1. Too high values do not provide enough structural information 
(see first rows in Figure 5). 



2. The appropriately informative value depends on the properties 
of the motion (thus of the microstructure) at a given location 
(compare columns of Figure 5, right side). 

3. Too low values force the isosurfaces to become extremely noisy 
(see last row of Figure 5). 

As the motion in highly organized tissue is less dispersed (i.e., 
a smaller support for P c f& which implies a larger support for S c fd 
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FIGURE 5 | On the left, one dimensional graphical representation of the 
choice of c-value. The drawing below the horizontal axis displays the 
structure of the sample in the infinitesimal volume element [— dx, dx] as 
seen by the magnetic moments from their initial position. The right side of 
the sample has less restrictive properties as on the boundary of tissue with a 
liquid, such as CC and cerebrospinal fluid (CSF). The noiseless isosurfaces 
consist of two points shown as dots on the graph. Low c-values correspond 



to the magnetic moments with longer travel paths providing more structural 
information than high c-values. However, too low values create noisy and 
disconnected isosurfaces represented with more than two points on the 
drawing. On the right, isosurfaces from different pixels in the baboon brain 
marked in Figure 3 demonstrate the effect of c-value on the information 
content from top = insufficient to bottom = noisy (see the text for the CSF 
column). 



thereby causing bigger truncation effects), increasing the diffu- 
sion gradient times 8 and A in section 2.3 will create almost "flat" 
displacement integral distributions at an isotropic medium like 
CSF. In this case, the small valued distribution [caused by con- 
stant integral value f -P c fd(X W d ) dW d = number of particles] is 
susceptible to noise, creating noisy isosurfaces of Figure 5 for all 
level values. In contrast, for the experiments conducted with an 
isotropic (water) phantom at much lower 8 and A values, the 
isosurfaces were spheres for a wide range of c-values (not shown). 

Figure 6 presents different isosurfaces that elucidate the tis- 
sue structure on the pons and the mid-brain of the fixed baboon 
brain sample described in section 2.4. The tracts on the left and 
right side of the mid-brain are visible with the ellipsoidal looking 
isosurfaces. Five isosurfaces from the same row of pixels marked 
on Figure 6 are displayed on Figure 6 corresponding to two dif- 
ferent c-values. The green isosurfaces with larger level values are 
smoother and less informative than the red ones with smaller c- 
values. Different viewpoints at each row of Figure 6 emphasize 
that the isosurfaces are 3-D objects. The figure demonstrates one 
of the challenges of presentation: the displacement integral val- 
ues must be considered in R 3 to grasp the complete information 
offered by CFD-MRI. 

Overall, the isosurfaces are not constrained to given forms like 
Gaussians, spherical harmonics or to any expansions. In fact, they 
are typically not even symmetric. They are structureless, general 
and direct. 

Isosurface visualizations constitute only one method to 
present the high dimensional information obtained from CFD- 
MRI. Another example is the dimension reduction by means 
of computing the so called orientation distribution func- 
tion (ODF) (Wedeen et al., 2005) obtained from radial inte- 
grals. However, for CFD-MRI the ODF raises the concern of 



inadequately presenting the "outskirts" of the P c fd because the 
dispersion of the outgoing rays shown in Figure 7 jeopardizes 
the inclusion of the values further away from the origin (see also 
Figure 6). 

New methods, additional elaborate schemes such as color cod- 
ing for representation of three dimensional functions aimed at 
displaying relevant microstructural information of CFD are left 
for future studies. 

4. DISCUSSION 

4.1. COMPARISON WITH EXISTING METHODS 

From a fundamental point of view, guided by the microstructure 
that surrounds them, the molecules are displaced due to ther- 
mal energy whether they are in the scanner or not. All existing 
DW-MR methods are designed with the same goal in mind: the 
reconstruction of the propagator 4 that describes the displacement 
of the magnetic moment from the DW-MR signal. 

However, as CFD-MRI demonstrates, from a systems science 
perspective the MRI scanner acts as a time-delay linear system 
with the input w, [displacement in Equation (4)], and the out- 
put Wi [displacement integral in Equation (8)], in Figure 8. The 
parameters A and 8 define the delay and filter parameters. Special 
attention is paid in CFD-MRI to isolate the Fourier variable, 
fc D = Gd from these parameters in contrast to the q-space vari- 
able (Callaghan, 1991): q = (27t) _1 Y8GD and the fo-value of 
DTI: b = y 2 II Gd II | & 2 (A - 8/3) (see Appendix A). 

The inverse problem of obtaining the propagator from the 
distribution of the displacement integrals is singular because of 



4 Also named nuclear spin self-correlation function (Callaghan et al., 1988) 
and translation probability (Cory and Garroway, 1990). 
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FIGURE 6 | The isosurfaces (P cfd = 0.14) on the pons and the 
mid-brain. Each of the boxes indicate an isosurface presented, 
respectively on the right. Each column presents the same isosurface 
from different view angles. The dots on top of the frames are placed for 



orientative purposes. The surfaces are not necessarily ellipsoids and they 
have mostly an asymmetric structure. The outer red surface is the set 
P cfd = 0.14 and the green surface is P cfd = 0.21. The red surface 
envelops the green one. 





FIGURE 7 | An example of a set of rays from the origin along 
which the distribution function is integrated to obtain orientation 
distribution function. As the rays disperse with increasing distance from 
the origin the points describing the displacement of a smaller number of 
particles contributes less and less to the numerical integration due to 



sparse sampling on both surfaces. The encapsulation of the isosurfaces 
on the right with larger level values by the smaller ones in Figure 6 
shows the information that would be missed with numerical radial 
integration. The utilization of isosurfaces is more informative as discussed 
in Figure 5 



the one-to-many relationship between the displacement and its 
integral: 

Wf = jw,(x)dx-j w,(x)dx Wi(f), (26) 

because all of the displacements with the same low frequency 
content in time are mapped to the same displacement integral 



value. This statistical accumulation prevents the determination of 
the propagator in a general environment from the distribution of 
displacement integrals 5 . 



5 One exception is the basic case of isotropic diffusion analyzed in Appendix A 
where the diffusion coefficient that describes the Gaussian propagator can be 
recovered using the adjustment factor for the displacement integral covariance 
namely b t in Equation (24) and the fo-value in DTI formulation. 
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FIGURE 8 | The DW-MRI signal from a signals and systems 
perspective (Ozcan, 2010b) describing the MRI observables of the 
diffusion phenomenon. The input is the magnetic moment displacement 
Wj and the output is the displacement integral, , defined in Equation (8). 



Existing methods' attempts to estimate the propagator relies 
on the narrow pulse approximation by assuming negligible pulse 
duration (8 = tdi — tdi = f</4 — in Figure 1) compared to 
pulse separation time, A, i.e., (8 << A) =>• (8 — > 0) specifi- 
cally in (A - 8/3) (Callaghan, 1991, p. 342). Under this short 
integration time assumption, in Wedeen et al. (2005) it is fur- 
ther argued that the approximation, wf ~ (xKte) — *i(fdi)) 8 = 
(Wj(to) — w,(£di))8, is plausible. Although by the intermediate 
value theorem and the sample path continuity of the Brownian 
motion (Shiryaev, 1995), the values of each integral in wf are 
attained at a time point within the respective integration intervals, 
[tdi , tdi\ and [te , tdi], the nowhere-differentiability of the sample 
paths (Shiryaev, 1995) implies that the intermediate time points 
satisfying the equality are not fixed as tdi and t#$, but are them- 
selves random variables. In consequence, without the inference of 
the displacements at fixed time points the propagator cannot be 
reconstructed. 

Moreover, elaborate derivations carried in Wedeen et al. 
(2005) to model the propagator as a conditional probability, 
p(x(td3)\x(tdi)), describing a Markovian process raise concerns 
specially in environments such as biological tissue since the 
particle's past collisions with microstructure guides its future dis- 
placements. In fact, while it violates the conditions of Wiener 
process (see Appendix A), this displacement memory provides the 
inference of the microstructure by way of affecting the displace- 
ments and consequently their displacement integral distributions. 
Accordingly, in CFD-MRI is indifferent to memory properties 
by modeling as a joint distribution function of random 
variables instead of the conditional probability of a stochastic 
process. 

A summary of CFD-MRI's detailed comparison with existing 
methods (Ozcan, 201 1; Ozcan et al., 2012) is presented below for 
completeness of exposure. Namely, there exits two avenues for the 
path from DW-NMR to DW-MRI in the literature: 



Mareci, 2003; Liu et al., 2004), and diffusional kurtosis imag- 
ing (DKI) (Jensen et al, 2005). 
2. Spectral methods originating from Callaghan's q- 
space (Callaghan et al., 1988) followed by the diffusion 
spectrum imaging (DSI) (Wedeen et al., 2005) and Q-ball 
imaging (Tuch, 2004). 

With the exception of the GDTI presented in Liu et al. (2004) 
[see also the discussion in Ozarslan et al. (2009)], all of the 
DW-MRI methods estimate symmetric quantities. The model 
matching methods project the data onto symmetric structures, 
such as ellipsoids in DTI or spherical harmonics in HARDI. The 
spectral methods use the magnitude of the signal in the Fourier 
transform resulting in symmetric functions (see section 2.2). It is 
difficult to imagine that molecular motion in a biological envi- 
ronment populated with different types of fluids, barriers and 
tissue would be symmetric at any given location, e.g., at the fiber 
junctions. Symmetry or lack of it ought to be determined by 
the data free of any constraints imposed by the model as in the 
implementation of CFD-MRI's unconstrained structure. 

The Fourier relationship between signal and joint distribution 
function provides a complete understanding of model matching 
methods. The methods start by applying DFT to the data in the 
first Z mr (imaging) dimensions. Thus, in the analysis of model 
matching methods the first Z mr independent variables are the phys- 
ical location. The three remaining untransformed variables are the 
independent variables of the Fourier reciprocal of displacement 
integral space, i.e., they are in the Fourier domain. The goal of 
model matching methods is to fit the preferred model to the dis- 
placement distribution function's Fourier transform, sampled at 
the (vector) values defined by the diffusion sensitizing gradients. 
In the case of Brownian motion this mixed variable (physical and 
frequency variables) approach is applicable because the diffusion 
coefficient D can be directly estimated from the Fourier domain 
by Equations (24) and (25). The mixed space, which works well 
for DW-NMR signal peak attenuation, is translated to DW-MRI 
at each position x by 



\ir^\ x ,k D )\ 

A-mr 



|/r mpleX (*,0)|exp(-H(fc D )) 



(27) 



where /™ m P lex j s gi ve n in Equation (18) and the function H 
defines the model, e.g., the quadratic form of DTI, spherical 
harmonics of HARDI or higher tensor expansions of GDTI 
[see Ozcan et al. (2012) for a detailed exposition]. In a sense, these 
methods' aim could be summarized as expanding the Fourier por- 
tion of the mixed signal space. The basic example with a single 
term in the expansion is DTI for which: 



H(ko) 



Y 2 b t klDk D 



(28) 



1. Model matching methods initiated by diffusion tensor imag- 
ing (DTI) (Basser et al., 1994; Mattiello et al, 1994) and 
further expanded with high angular resolution DW imaging 
(HARDI) (Frank, 2001), composite hindered and restricted 
model of diffusion (CHARMED) (Assaf and Basser, 2005), dif- 
fusion orientation transform (DOT) (Ozarslan et al., 2006), 
two versions of the generalized DTI (GDTI) (Ozarslan and 



where the calculation of b t from the PDE approach in Ozcan 
(2010a) and covariance of displacement integrals in Appendix A 
resulted in the same value: S 2 (A - 8/3). In CFD parlance, by 
the Fourier relationship between Gaussians in Equation (25), the 
diffusion quadratic form, D, is estimated in /cd -space, without 
recourse to a Fourier transform because of its direct appearance 
in Equation (28) rather than its inverse, D _1 , in the Gaussian 
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of motion space. The coefficient matrix defined by carefully 
selected vectors in -space that satisfy the invertibility condi- 
tions (Ozcan, 2005) is used for the estimation of D in the linear 
algebraic framework of symmetric matrices (Papadakis et al., 
1999; Ozcan, 2010a) [also refer to Ozcan (2010a) and Ozcan et al. 
(2012) for the correspondence with the B-matrix formulation 
ofBasser et al. (1994)]. 

The magnitude-based Fourier relationship presented in q- 
space methodology (Callaghan et al., 1988) is the origin of 
spectral methods. In Callaghan's book (Callaghan, 1991), par- 
allel to the historical development of DW models, the the- 
ory is first developed for NMR experiments [see Callaghan 
(1991, Chap. 6)] using polarized neutron scattering analogy. 
However, the translation from NMR to MRI is presented 
[see Callaghan (1991, Chap. 8)] asserting without proof that the 
imaging and displacement portions of the signal are separable 
[see Callaghan (1991, Chap. 8, pp. 440)]. The derivations of 
section 2.1 demonstrate that this is not the case. 

In addition, by the affine dependence of /c rw on k mr CFD 
derivations show that the inseparability partially accounts for the 
non-Hermitian nature of the S c fd- Taking the magnitude of the 
DW-MRI signal, as in the case of DSI (Wedeen et al., 2005), does 
not count as a phase correction. Figure 9 demonstrates that by 
preserving Hermitian property, CFD-MRI captures correctly the 
crossing fibers at the junction of the CC and EC. 

4.2. CONCLUSION AND FUTURE STUDIES 

In the biomedical imaging modalities' grand aim of biomarker 
capability establishment, the discovery path for CFD-MRI passes 
through the distribution function: 

5 c fd -> Pdd -> Biologicalproperties. 

With P c ^ in the middle, both sides of the path present themselves 
with important challenges. 

First and foremost, in DW-MRI, the displacements with- 
out reference to initial positions [see Equation (45)] prevent 
the inference of microstructure position. For example, the 



distribution function of the biological phantom (Ozcan et al., 
2011) constructed with two crossing rat trigeminal nerve fibers is 
always in the form of two crossing bars across the origin regard- 
less of the nerves' position as long as their relative angle is kept 
the same. Also in the same phantom, the agar gel (isotropic com- 
ponent) appears as a sphere around the origin of domain 
without the possibility of identifying its location. As the distribu- 
tions from various types of microstructural components accumu- 
late around the origin, the discrimination level of overlaps, more 
prominent with increasing biological tissue complexity, directly 
defines the sensitivity and specificity for microstructure changing 
pathologies. The important goal is the assessment of the theoret- 
ical aspects of the distribution function in order to understand 
whether it can detect in a timely manner, e.g., before signifi- 
cant disease progression, those changes. The determination of 
biophysical conditions behind the asymmetry (see also Ozarslan 
et al, 2008; Ozarslan, 2009) in the distribution functions is also 
part of the same goal. 

However, the absence of analytical descriptions for P c fd even 
in simple environments requires the investigations to be con- 
ducted with numerical simulations (Ozcan et al., 2011) of par- 
ticle motion within carefully designed geometries (Landman 
et al, 2010) and locally variable diffusivities. Along with numer- 
ical phantoms mimicking biological ones (Ozcan et al., 2011), 
histopathological information is also being used for interpreta- 
tion and validation (Budde and Frank, 2012; Budde and Annese, 
2013). Additionally varying the time parameters 8 and A will 
exploit the filtering effects caused by the displacement integral 
that will determine whether further information extraction is 
possible by expanding data acquisition with an appropriate set of 
parameter values. 

On the other hand, on the discovery path's initiation by CFD- 
MRI signal formation, the re-establishment of Hermitian sym- 
metry requires, in addition to the theoretical reasons presented 
herein, the analysis and quantification of Hermitian disruptive 
artifacts and systematic conditions in real data. Constructed by 
initially experimenting with elementary phantoms (e.g., water 
and mineral oil), this signal model expansion is necessary for 



CC Junction White Matter Corpus Callosum CSF 




FIGURE 9 | The comparison of P cfd (Ozcan, 2011) (first row) with the 
diffusion spectrum imaging orientation distribution function (DSI-ODF) 
(second row). Both functions are calculated from the same data on the 
right junction of the corpus callosum (CC) and the external capsule (EC), 



specifically from the pixels marked on the Figure 3. The isosurface 
(P cfd = 0.17) captures the asymmetric structure of the fiber crossings while 
the ODF is constrained to be symmetric for all of the pixels. Note that in CSF, 
ODF detects structure which is not present in reality as indicated by CFD. 
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the development of more elaborate systematic phase corrections, 
possibly by utilizing complex analysis theory. Specifically, accu- 
rate estimation of the pertinent Fourier transform of P c fd from 
real data points in a clinical setup is targeted by the adaptation 
of CFD-MRI to fast sequences 6 , such as echo planar imag- 
ing (EPI) prone to Eddy current artifacts. The model will be 
expanded up to the point of reaching only minimal incremental 
improvements with new phase correction algorithms. Thereafter, 
relying on the residuals' content, which is free from displacement 
effects consequent to the application of system-wide uniform 
phase corrections 7 , more effective biomarker construction would 



Acquisition time is shortened with reduced sampling schemes aimed at spe- 
cific goals, e.g., compressed sensing for tractography (Landman et al., 2012) 
at the expense of some information loss. 

7 Instead of pixel by pixel corrections that would completely eliminate the 
imaginary part in Figure 3. 



be possible by the inclusion of extra information such as tis- 
sue susceptibility (Liu et al., 2013). Likewise, on a larger scale 
CFD-MRTs general aim is to improve outcomes of multimodal 
imaging, e.g., in prostate cancer strategies (Turkbey and Choyke, 
2012). 
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GLOSSARY 

w;: ith magnetic moment's displacement from its initial position 

Xi{t 0 ). 

W 1 : The difference of magnetic moment displacement inte- 
grals during the periods diffusion (d) gradients are turned 
on. Its distribution is the main quantity of interest provid- 
ing information about micro structure. Joined with the other 
displacement integrals during acquisition and rewind times, 
W acq , W rw , respectively, it forms the total displacement integral: 
W = (W d , W™, Wf) € R 7 . 

J'cfdC*) W 1 ): CFD joint distribution function of the number of 
magnetic moments with initial position x possessing displace- 
ment integral values of W d . It is obtained by marginalizing 



□total 
cfd 



(*, W). 



ScfdCfcmr, fco, fcrw): DW-MRI signal that comes out of the scanner 



which is equal to the Fourier transform of P?S . Theoretically it 
must be Hermitian symmetric since P]^ 1 is real valued. 
(k mT , kj)): fc c fy-space variables corresponding to imaging gradi- 
ents (mr) and diffusion (D) gradients, respectively, 
in,: Rewinding (rw) frequency variable affine dependent on k ml . 
It cannot be sampled at 0 because of its dependence. For P^" 1 
marginalization, phase corrections are applied to estimate S c fd at 
/c rw = 0. 

jcompkx^ Complex valued images obtained by taking the 
Fourier transform of S c a with respect to imaging frequency k ml . 



Ikm*( x > ^d)'- Image domain phase corrected t 



complex 



(x, k D ). 



APPENDICES 

A. SECOND ORDER STATISTICS FOR DISPLACEMENT INTEGRALS OF 
SELF DIFFUSION 

In this section, the covariance of displacement integrals is 
derived for the special case of particles executing Brownian 
motion. The following preliminaries are necessary to calcu- 
late the second order statistics for the displacement integrals 
of Equation (9). It is assumed for ease of notation below that 
tm — t m + l ■ 

The Wiener process (Shiryaev, 1995), w, which describes 
the diffusion in an isotropic homogenous medium, is sam- 
ple path continuous, has w(0) = 0 and independent incre- 
ments (i.e., Markovian property) with a normal distribu- 
tion, meaning w(t) — w(s) ~ N(0, (f — s) D) where t > s > 0 
and D > 0. Using these properties, the covariance of the 
displacement is: 

E [w(f) w r (s)] = E [(w(f) - w(s) + w(s)) w r (s)] 

= E \w(s)w T (s)\ = min(f, 5) D. (29) 

It is straightforward to show that the mean for the displacement 
integral processes is 0. The covariance in the case of a single 
interval is calculated as 



' 12 

/ 



W{(x) dx 




wj(s)ds\ 



H t 2 



j j e[w,(x) wj '(5) J dxds 



h h 

// 



1 7 

min(T, s)Ddxds = -(t 2 - ti) (2fi + t 2 )D, (30) 



and for a non-overlapping interval: 



j j E j^Wj(x) wf(s) j dxds = j j min(T, s) D dx ds 



f3 fl 
1 



f 3 ti 



(t 4 -t 3 )(4-tj)D. 



(31) 



Finally, using the formulas above, the time points of Figure 1, 
[ti ti t 3 = [tji tfo tdi t^], and these substitutions, 

8 = t 2 — fl = frf2 — fdl = 1 4 — f 3 = td4 — t<te, 
A = t 3 — fi = fj 4 — frf 2 = f 4 — f 2 = te — f<il . 

the following is obtained for Wf : 



Wf 



(2f 3 -r-f 4 ))/3-(f4-f3)(f 2 2 -ff)]l> 



8 2 (A-8/3)D. 



(32) 



The scalar factor, 8 2 (A — 8/3), defined as fo t -value in Ozcan 
(2009, 2010a), is a function of time directly related to the 
measured quantity, the displacement integrals. It is completely 
detached from the measurement parameters such as diffusion 
gradient strength, in contrast to the fo-value used in the lit- 
erature, b = y 2 || Gd II 1 b t . The dependence of the derivations 
on Markovian property restricts the validity of the calcula- 
tions to isotropic samples where the diffusion is character- 
ized by the diffusion coefficient, D, in the distribution of 
magnetic moment displacement w. For an isotropic homoge- 
nous medium, D is estimated from the displacement inte- 
gral's (W d ) covariance b t D= 8 2 (A - 8/3) D using DW-MRI 
acquisitions. 

The rigorous treatment of the theory of the stochastic pro- 
cesses in Doob (1990, p. 62) demonstrates that the sample paths 
of stochastic processes are Lebesgue integrable under certain 
conditions and these integrals are well defined random vari- 
ables. Although a rigorous mathematical analysis, which proves 
that these conditions are met, is out of the scope of this 
manuscript, it is important to note that Equation (32) does 
not prove that the displacement integrals are Gaussian ran- 
dom variables. In this case, the central limit theorem does not 
apply because the displacement integrals of Equations (8, 9), 
are not the sums of independent variables since in the integral 
approximating sum 

rh N 
I Wj{x) dx ~ lim Y w,(fo + k dx) dx 
Jt 0 N ^°°k=l 



(N 00, dx -+ 0, N dx = ti - f 0 ), 



(33) 
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Wi(x), rather than independent increments, Wi(t) — Wi(s), is 
present. A variant of central limit theorem for sums of dependent 
variables satisfying certain conditions (Shiryaev, 1995, p. 541) 
might provide the theoretical validation of this highly theoretical 
open problem. 

B. DERIVATION OF THE PHASE EQUATION 
B1. Bloch equations and their solutions 

The time evolution of the magnetization, m,(f) e R 3 , that gener- 
ates the MR signal is modeled by the Bloch equations: 



by calculating the matrix exponential using f2(x,(r), ro) 
f* o oi(xi(x))dx = f* to G{Xi, x) ■ Xi{x)dx: 



J to 



exply / A(x)dx 



cos(yf2 (x,(t), t 0 )) sin(yfi (x,(t), t Q )) 0' 
-sin(yC2(x,(f), to)) cos(yQ(Xi(t), t 0 )) 0 
0 0 1 



Equation (38) defines a system with a rotating transverse magne- 
tization component that can be written in a complex form w; = 
m ix + J m iy The corresponding first order differential equation is 



dm,(t) 
dt 



y ntiit) x B(xj, t) 



l 

r 2i 



0 0 
i 



0^-0 

0 0 TfJ- 



m,(f) 



+ 



1 



0 
0 



(34) 



with Tu and T% denoting spin-lattice and spin-spin relaxation 
constants, respectively. In the MR scanner the magnetic field as a 
function of time and position is given by 



(35) 



where Bo is the static (strong) magnetic field, B" ! (f) is the radio 
frequency (RF) pulse modulated at the precession frequency of 
Bo and G(x, t) e ]R 3 describes the magnetic field gradients. The 
expression of the magnetic field simplifies further in the rotating 
frame: 





' 0 " 




0 


B(x, t) = 


0 


+ Bf(0 + 


0 




_B 0 _ 




. G(x, t) 



B(x, 0=Bi(t) + 



0 
0 

G(x, t) ■ x(t) _ 



(36) 



When the magnetic field gradients are turned on without the RF 
pulse the system matrix becomes 





0 




0 


u>(xi(t)) 0" 




fflj X 


0 




-oo(x,'(f)) 


0 0 


m, = A(t) mi. (37) 




_w(x,(f))_ 




0 


0 0_ 





with oo(Xj(£)) = G(xi, t) ■ Xi(t). As the time-dependent system 
matrix A(t) commutes with its integral f A(x)dx, when the 
relaxation terms are negligible the solution of Equation (34) is 
obtained [see Rugh (1995, pg. 59)] as: 



m;(f) = exp 



A(x)dx m,(fo) 



(38) 



dnii(t) 
dt 



-JY<a(.Xi(t))mi(t), 



with the solution 



m ; (f) = exp(-j yQ(Xi(t), t 0 )) m,(fo)- 



(39) 



(40) 



For the portion of the pulse sequence involving the RF pulses, 
in general (Bernstein et al, 2004) the magnetic moments are 
assumed to be immobile [i.e., Wi(f) = 0 in Equation (4)] thereby 
defining a constant system matrix that has a standard matrix 
exponential that solves the differential equation. Herein, however, 
the effect of the displacements specifically during the it -pulse 
results in a time varying system matrix. Without loss of gen- 
erality, the jt -pulse is applied in the direction 8 [1, 0, 0] as 
in Hinshaw and Lent (1983), along with the slice select gradi- 
ent, Gn = [0, 0, ^^3], active only during the RF pulse. The RF 
pulse modulated at the same frequency corresponding to the cen- 
ter of the slice, Bo + g U 3 x c i, becomes a constant vector, Bi(f) = 
[Bi a , 0, 0] J , in Equation (36). The system of equations in this 
rotating frame is 



ntj x 



Bla 

0 

a>jt(*;(0). 



m,- = A(t) mi 



0 G>n(Xi(t)) 0 

-<A>„(Xi(t)) 0 B la 
0 -B lfl 0 _ 

(41) 

with 9 mJMt)) = G„ • (xi(to) + Wi{t) - x c ) = g n3 {x ii (t 0 )+ 
W3i(0 —x c i). The procedure would be to solve Equation (41) 
and then express the solution in the initial rotating frame to 
obtain a phase component, $ m -, that would be a function of 
the displacement integral and the initial position from the slice 
center. 

However, the new system matrix, A(t), does not commute with 
its integral preventing the calculation of a matrix exponential for- 
mulating a convenient analytical solution as in Equation (40). 
With the rigorous analysis describing the exact manifestation of 
the RF pulse left to future studies, herein due to the large magni- 
tude of the RF field in comparison to the magnetic field gradient, 
its effect will be approximated as a sign reversal of the phase. 



8 The choice of direction does not matter (Haacke et al., 1999, p. 387). 
9 The absence of a rewinding slice select gradient, in contrast to it/2-pulse 
in Figure 1, requires the addition of the slice center, x c , in the formulation 
as there will be a phase shift along the slice select direction in the slab [see 
Equation (11) and footnote 1 ] . 
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In this work, the phase change induced by the magnetic field gra- 
dient is considered as a factor that is automatically corrected by 
the phase correction algorithm of section 2.2. 

On the other hand, the effect of the displacement induced 
by the RF slice select gradient (SS-EX in Figure 1) during the 
excitation period is systematic. As both the RF pulse and the mag- 
netic field gradient remain fixed, statistically speaking, the phase 
of the total transverse magnetization is affected by the displace- 
ments in the same amount at each acquisition prior to the tipping. 
The position and displacement encoding occurs thereafter dur- 
ing the remainder of the pulse sequence (see also footnote 1 for 
the phase correction along the slice select direction within the 
slab) leaving the jt/2-pulse out of the derivation of the phase 
equations below. 

B2. Phase equation 

For clarity of exposition, the calculations are carried out with 
the assumption of ideal gradient amplifiers creating rectangular 
shaped gradient pulses in the time domain and perfect linear- 
ity in physical space. In short, during the time interval when the 
amplifiers are powered on G(x, f) = G € R 3 . 

The derivation of Equations (5- 7) provided in Ozcan (2012) 
is briefly presented in this section for reference purposes. Using 
the definitions of the variables in Figure 1, the evolution of the 
phase is described as: 

1. First, imaging gradients for read out rewinding, G ro e M 3 , 



phase encoding, G p 



and slice select G ss € R on the 



time interval [to, f rw ] are turned on. After these gradients are 
applied the phase and the magnetization become 

* rw 

= j (GroW + Gp e (T) + G ss (t)) • Xi(x) dx 
to 

( 1 \ 

= (G ro + Gpe + G ss ) • Af rw x,(f 0 ) + / Wi(x) dx I 



Af lw (G ro + Gpe + G ss ) ' Xi(t Q ) 

'rw 

+ (G ro + Gpe + G ss ) • J w,(t) dx. 
OTi(frw) = exp(-; y ^™) nii(t 0 ). 



(42) 



(43) 



2. The RF jt-pulse between the diffusion gradient pulses, Gd G 
K 3 , and the accompanying slice select gradient provide the- 
oretical sign reversal of the phase. The slice select gradient's 
encoding of magnetic moment motion into the signal is 
expressed by <t> n ; (see Appendix Bl). Since = x;(fo) + 

Wi(tdk), k= 1, ...,4; atf = t d4 



+ j G D -x,(x)dx 



1,12 



- j G D -x,(x)dx (44) 



<t>m + G D • (frf4 - td3)xi(t Q ) + / Wi(x) dx 



tdi 

o) + /v 



\ tdi / 

- G D • ^(t d2 - fdi)x,(f 0 ) + j w,(t) ckj 

(tdi tj 2 \ 

j wi(x) dx- j w,{x)dx\ 
di t dl } 

+ ((tdi ~ te) - (tdi ~ tdi)) G D • Xi(fo). (45) 

If the diffusion gradient times are equal, i.e., td4 — 
tdj, = tdi — tdi = 8) then the last term in Equation (45) 
is equal to zero, erasing the influence of initial posi- 
tion from motion encoding part of the signal. This is 
the insight gained by using the formulation with dis- 
placement integrals, f Wi(x) dx, rather than the center of 
mass (COM) of random walk, J x(x) dx, introduced in 
Mitra and Halperin (1995). 

The sign reversal affects also previously accumulated phase: 

m(tdi) = exp(-j y ^d) exp(; ^ rw ) m(to)- (46) 
3. The last part of the sequence is where the data are collected: 

^ro(f) = / G ro • (x,(t 0 ) + w,(t)) dx = (t - f acq ) G ro 



/ 

'acq ( 
■ Xi(t 0 ) + j 



G ro • wi(x)dx, 



(47) 



leading to 

m,(f) = exp(-; y (f2i-o(0 + fi D - ^rw)) W;(f 0 ). (48) 
and Equations (5- 7). 

C. TOTAL NUMBER OF PARTICLES FROM CFD-MRI 

The fundamental Fourier relationship of Equation (16), 
Scfd(kmr, fc D, U = ^{P^f l }(k m i, k D , /c rw ) establishes the rela- 
tionship between the standard MR image space and the higher 
dimensional CFD space. By Equation (16) 

Scfd(fcmr, 0, 0) = J P%f l (x', W) exp [-J Yfcmr ■ *'} dx' dW 

(49) 

and by Equation (18), f k ° mplex (x, k D , k, w ) = T7 l {S c fd(k mr , 
)}, the image obtained without the diffusion encoding 
magnetic field gradients is 



7 complex ( ^ 0> Q) = j T-l {Scfd(femrj 0> Q)} 



(50) 



/ 



P^(x',W) exp{-j Y fcn 



• (x' - x)} dk mr dx' dW 
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by the Fourier Property : 



/ 



exp{± j 2tt kx] dk = S(. 



(51) 



~ total number of particles at x. 



dW 



(52) 

(53) 
(54) 
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